clear all
%% Chosing stocks in the same sector and equivalent liquidity
%Ex: in Industry sector(sector_id = 2)
all_sec_id = get_info_from_database('security-sector-2', 20);

idx_ok = [];
med_turnover_pct = [];

for i = 1:size(all_sec_id,1)
    rslt = get_info_from_database('daily_stats', 'sec_id', all_sec_id(i));
    if ~isempty(rslt)
        idx_ok = horzcat(idx_ok, i);
        med_turnover_pct = horzcat(med_turnover_pct, rslt.med_turnover_pct);
    end
end

[vals idx] = sort(med_turnover_pct,'descend');
all_sec_id_ok = all_sec_id(idx);
%% Get financial ratios from selected stocks
sub_sec_id_ok = all_sec_id_ok(:,1);
financial_ratio = [];
idx_ok = [];
for i = 1 : size(sub_sec_id_ok,1)
    rslt = get_info_from_database('financial_ratio', 'sec_id',...
        sub_sec_id_ok(i),'from', 2008,'to',2011);
    if ~(iscell(rslt) && strcmpi(rslt{1}, 'No Data'))
        idx_ok = horzcat(idx_ok, i);
        financial_ratio = vertcat(financial_ratio, rslt);
    end
end
sec_id_ok = sub_sec_id_ok(idx_ok);
financial_ratio = financial_ratio(:,[1:4 16:20 24 26]);
financial_ratio = cell2mat(financial_ratio);
dataX = [];

for j = 1:size(financial_ratio,1)
    if sum(isnan(financial_ratio(j,:)))== 0
        dataX = vertcat(dataX, financial_ratio(j,:));
    end
end
%% Normalize data X
% Compute the quantile
[value1 idx1] = sort(dataX(:,5:end));
parameters = dataX(:,5:end);
all_q_X = [];
for i = 1 : size(parameters,2)
    q_X(idx1(:,i)) = (1:length(parameters))/length(parameters);
    all_q_X = horzcat(all_q_X, q_X');
end
dataX_normalize = horzcat(dataX(:,1:4),all_q_X);

%%
mat_X = [];
vec_y = [];
datequarter=[];
t1=1;t2=1;
res=[];IndexPos=1;firstprice=[];
for y = 0:3
    for q = 0:3
        year = 2008+y;
        quarter = 1+q;
        indLoc = intersect(find(dataX_normalize(:,2)==year),...
            find(dataX_normalize(:,3)==quarter));
        if ~isempty(indLoc)
            sec_id = dataX_normalize(indLoc,4);
            ret = [];
            if quarter ==1
                fromTime = sprintf('30/06/%d', year);
                toTime = sprintf('30/09/%d', year);
            elseif (quarter==2)
                fromTime = sprintf('30/09/%d',year);
                toTime = sprintf('31/12/%d',year);
            elseif (quarter==3)
                fromTime = sprintf('31/12/%d',year);
                toTime = sprintf('31/3/%d',year+1);
            else
                fromTime = sprintf('31/03/%d',year+1);
                toTime = sprintf('30/06/%d',year+1);
            end
            for i = 1:length(sec_id)
                result = get_info_from_database('security', 'sec_id', sec_id(i) , ...
                    'from', fromTime, 'to', toTime);
                if isempty(result)||length(union([],result(:,6)))/length(result(:,6))<0.5
                    indLoc(dataX_normalize(indLoc,4)==sec_id(i))=[];
                    continue
                else
                    result1 = (result(2:end,[2 4 5 6])-result(1:end-1,[2 4 5 6]))./result(1:end-1,[2 4 5 6])*100;
                    ret=[ret;mean(result1(:,4))];
                    datequarter=[datequarter;size(result1,1)];
                    if isempty(res)||~ismember(sec_id(i),firstprice(:,5))
                        firstprice(t1,1:4)=result(1,[2 4 5 6]);
                        firstprice(t1,5)=sec_id(i);
                        t1=t1+1;
                    end
                    ttt=size(result1,1);
                    res(IndexPos:IndexPos+ttt-1,1:4)=result1;
                    res(IndexPos:IndexPos+ttt-1,5)=repmat(sec_id(i),ttt,1);
                    IndexPos=IndexPos+ttt;
                    result1=[];
                    result=[];
                end
            end
            dataXyq = dataX_normalize(indLoc,5:end);
            mat_X = vertcat(mat_X, dataXyq);
            vec_y = vertcat(vec_y, ret);
            
        end
    end
end
[res dd]=sortrows(res,5);
[pos,~]=secPos(res(:,5));
firstprice=sortrows(firstprice,5);
[b,~,~,inmodel] = stepwisefit(mat_X,vec_y);
b=b';
YFIT = b(inmodel)*mat_X(:,inmodel)';
YFIT=maniArray(YFIT,datequarter);
YFIT=YFIT(dd);
mat_XX=[];
window=0;zz=[];
vec_x=[];
zzz=[];
for i=1:length(firstprice(:,2));
    i
    price=[];
    
    retvet=res(res(:,5)==firstprice(i,5),1:4);
    volume=CalPrice(firstprice(i,1),retvet(:,1));
    price(:,3)=CalPrice(firstprice(i,4),retvet(:,4));
    price(:,1)=CalPrice(firstprice(i,2),retvet(:,2));
    price(:,2)=CalPrice(firstprice(i,3),retvet(:,3));
    price = adjPrice(price);
    [mat_X window] = patternRec3(price(2:end-1,3));
    mat_XX=[mat_X price(1+window:end-1,3)];
    
    % ADX indicator
    [vout]=indicators(price(2:end-1,:),'adx',3);
    mat_XX = [mat_XX vout(window:end,3)];
    [vout]=indicators(price(2:end-1,:),'adx',14);
    mat_XX = [mat_XX vout(window:end,3)];
    vout=[];
    
    % Stochastic Oscillator (Base on
    [vout1]=indicators(price(2:end-1,:),'fsto',3,3);
    mat_XX = [mat_XX vout1(window:end,:)];
    
    [vout2]=indicators(price(2:end-1,:),'fsto',14,3);
    mat_XX = [mat_XX vout2(window:end,:)];
    mat_XX = [mat_XX vout1(window:end,:)./vout2(window:end,:)];
    vout1=[];vout2=[];
    
    % momentum (rate of change)
    [vout1]=indicators(price(2:end-1,3),'roc',3);
    mat_XX = [mat_XX vout1(window:end)];
    
    [vout2]=indicators(price(2:end-1,3),'roc',14);
    mat_XX = [mat_XX vout2(window:end)];
    
    mat_XX = [mat_XX vout1(window:end)./vout2(window:end)];
    vout1=[];vout2=[];
    
    % relative strength index
    
    [vout1]=indicators(price(2:end-1,3),'rsi',3);
    mat_XX = [mat_XX vout1(window:end)];
    
    [vout2]=indicators(price(2:end-1,3),'rsi',14);
    mat_XX = [mat_XX vout2(window:end)];
    
    mat_XX = [mat_XX vout1(window:end)./vout2(window:end)];
    vout1=[];vout2=[];
    
    % macd
    [vout1]=indicators(price(2:end-1,3),'macd');
    mat_XX = [mat_XX vout1(window:end,:)];
    
    %average true range
    [vout1]=indicators(price(2:end-1,:),'atr',3);
    mat_XX = [mat_XX vout1(window:end)];
    
    [vout2]=indicators(price(2:end-1,:),'atr',14);
    mat_XX = [mat_XX vout2(window:end)];
    mat_XX = [mat_XX vout1(window:end)./vout2(window:end)];
    vout1=[];vout2=[];
    
    % simple moving average
    
    [vout1]=indicators(price(2:end-1,3),'sma',3);
    mat_XX = [mat_XX vout1(window:end)];
    
    [vout2]=indicators(price(2:end-1,3),'sma',14);
    mat_XX = [mat_XX vout2(window:end)];
    mat_XX = [mat_XX vout1(window:end)./vout2(window:end)];
    vout1=[];vout2=[];
    % Calculate real residual
    ttt=length(retvet(:,1));
    if window>ttt   %This is to prevent the case of revet vector length lower than the window
        window=ttt;
    end
    if i==length(firstprice(:,2))
        YFIT(pos(i)+window:end)=YFIT(pos(i)+window:end)-retvet(window+1:end,4);
%         zz(i,2)=length(YFIT(pos(i)+window:end));
    else
        YFIT(pos(i)+window:pos(i+1)-1)=YFIT(pos(i)+window:pos(i+1)-1)-retvet(window+1:end,4);
%         zz(i,2)=length(pos(i)+window:pos(i+1)-1);
    end
    % Remove window part
%      qua dai dong doan nay/
    if size(mat_XX,2)==34;
        vec_x=[vec_x;mat_XX];
        YFIT(pos(i):pos(i)+window-1)=[];
    elseif i==length(firstprice(:,2))
        YFIT(pos(i):end)=[];
        window=window+length(YFIT(pos(i)+window:end));
    else
        YFIT(pos(i):pos(i+1)-1)=[];
        window=window+length(YFIT(pos(i)+window:pos(i+1)-1));
%         zz(i,1)=size(vec_x,1);-+
    end
        mat_XX = [];
        pos=pos-window;
    end